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Investigations  of  relevant  issues  involved  in  developing  a  reduced  order  model  (ROM) 
are  performed  for  describing  combustion  response  to  specific  excitations  using  Euler 
equations  as  a  continued  work  from  a  previous  studies  using  a  reaction-advection  scalar 
equation.  The  ROM  is  obtained  by  employing  Galerkin’s  method  to  reduce  the  high-order 
PDEs  to  a  lower-order  ODE  system  by  means  of  POD  eigen-bases.  For  purposes  of  this 
study,  a  linearized  version  of  the  Euler  equations  is  employed.  The  knowledge  obtained  from 
previous  scalar  equation  studies  is  applied  to  ROM  development  of  Euler  equations. 
Numerical  stability  issues  are  encountered  in  Euler  equation  studies,  the  cause  of  which  is 
narrowed  down  to  normalization  methods  for  vector  variables.  The  effects  of  normalization 
methods  are  then  further  assessed  in  terms  of  the  ROM  characteristics  and  its  predictive 
capability. 


I.  Introduction 

^jombustion  instability  is  a  complex  phenomenon  that  results  from  the  coupling  between  the  modes  of  heat  release 
and  acoustics.  In  practical  combustor  devices  the  complexity  is  greatly  amplified  by  turbulent,  compressible 
flow,  very  high  rates  of  heat  release,  and  often  complicated  geometries  and  acoustic  boundary  conditions. 
Modem  computational  capability  offers  the  potential  for  moving  beyond  the  empirically-based  design  analysis  of 
the  past,  but  simulations  of  full  scale  dynamics  for  engineering  analysis  are  still  out  of  reach.  However,  high  fidelity 
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simulations  of  smaller  scale  domains  can  be  used  to  obtain  reduced  order  models  of  the  combustion  response  that 
can  presumably  provide  accurate  descriptions  of  the  linear/nonlinear  coupling  between  acoustics  and  combustion. 

The  model  reduction  techniques  have  been  well  studied  to  develop  numerically  stable  and  robust  reduced  order 
models  (ROM)  [1-3].  ROM  techniques  have  wide  application  in  non-reacting  flow  problems  including  flow  control 
[4-6]  and  unsteady  aeroelasticity  [7,  8].  Recently  application  has  been  extended  to  combustion  problems  [9,  10].  A 
preliminary  exploration  of  the  POD/Galerkin  technique  using  a  model  reaction-advection  scalar  equation  in 
developing  valid  and  capable  reduced-order  model  was  performed  by  Huang  et  al.  [11],  which  mainly  assesses  the 
capability  of  the  ROM  for  predicting  responses  at  target  frequencies.  Other  than  the  implementation  of  techniques  in 
practical  problems,  focuses  have  also  been  placed  on  resolving  the  robustness  and  stability  issues  in  developing  a 
valid  reduced  order  model.  As  reported,  the  issues  can  come  from  the  possible  inherent  lack  of  numerical  stability  in 
the  POD/Galerkin  method  itself  [12],  truncation  of  low-energy  dissipative  POD  modes  [1]  and  simplifications  of 
higher-order  equations  [13].  The  balanced  POD  technique  has  been  proposed  to  build  numerically  stable  ROM  for 
linear  systems  [2,  14].  Bergmann  et  al.  proposed  to  add  residuals  of  the  Navier-Stokes  equations  to  account  for  the 
absence  of  low-energy  dissipative  POD  modes  [1].  Barone  et  al  [15,  16]  proposed  to  stabilize  the  reduced  system  by 
symmetrizing  the  higher-order  PDE  with  a  preconditioning  matrix.  Rowley  et  al.  also  pointed  out  that  defining  a 
proper  inner  product  can  be  important  when  dealing  with  model  reduction  of  the  Navier-Stokes  equations  [3]. 

In  this  paper,  we  extend  our  previous  exploratory  work  using  a  representative  scalar  equation  [11]  to  vector 
equations.  The  remainder  of  the  paper  is  organized  as  follows.  In  Section  II,  a  summary  of  previous  work  and  major 
conclusions  is  given.  In  Section  III,  we  present  the  Euler  equations  with  a  modeled  combustion  source  term  used  in 
the  present  study  as  well  as  the  Galerkin  formulation  and  POD  techniques  for  deriving  reduced  order  models  for  the 
linearized  version  of  the  equation.  In  Section  IV,  we  describe  the  test  problem  and  the  representative  characteristics 
of  the  CFD  solutions  used  in  the  ROM  studies.  In  Section  V,  we  follow  the  procedures  developed  in  our  previous 
scalar  equation  studies.  First,  we  validate  the  resulting  POD  eigen-bases  and  ROMs.  Second  a  specific  numerical 
stability  issue  is  identified  during  ROM  development  that  can  be  attributed  to  the  normalization  methods  used  in 
generating  the  POD  eigen-bases  for  vector  variables.  Based  on  that,  the  effects  of  normalization  methods  on  ROM 
development  are  further  evaluated.  In  the  final  section,  we  provide  some  concluding  remarks  and  future  directions 
for  continued  research  and  development. 


II.  Review  of  Scalar  Equation  Studies 

A  review  of  previous  work  using  a  reaction-advection  scalar  equation  [11]  is  given  in  this  section.  A  Galerkin 
approach  based  upon  Proper  Orthogonal  Decomposition  (POD)  techniques  is  used  to  assess  the  performance  and 
capabilities  of  the  derived  reduced  solver.  The  model  equation  is  posed  in  a  non-linear  form  but  linearized  versions 
are  used  also  to  allow  more  detailed  investigation  of  the  impact  of  source  term, 
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Model  Equation:  +  a6V’r)  =Kl\\^9(X,T)\6{X,T)  +  Kl0T\e(X,T)-0(X)\ 

dO'(X,z )  dd'(X,T 

Linearized  Version:  - - - -  + - - - 

dr  dX 

The  approach  uses  numerical  solutions  of  the  model  equation,  obtained  by  perturbing  quantities  of  interest  such 
as  the  inlet  conditions,  as  the  database  for  building  a  set  of  POD  eigen-bases.  The  POD-derived  eigen-bases  are,  in 
turn,  used  to  obtain  an  ordinary  differential  equation,  which  constitutes  the  ROM.  Once  the  ROM  is  established,  it 
can  be  used  as  lower-order  test-bed  to  predict  detailed  results  within  certain  parametric  ranges  at  a  fraction  of  the 
cost  of  the  full  governing  equations.  Major  conclusions  from  the  studies  using  the  linearized  model  scalar  equation 
are  summarized  below  and  further  applied  in  the  Euler  equations  studies  given  in  the  next  section: 

1 .  In  POD  eigen-bases  generation,  for  cases  where  the  numerical  rank  of  POD  matrix  is  equal  to  the  analytical 
rank,  the  number  of  POD  modes  is  twice  the  number  of  frequencies  included  in  the  forcing  function; 

2.  Consistent  discretization  between  the  CFD  solution  and  the  ROM  is  required  to  obtain  good  accuracy; 

3.  Including  the  same  number  of  POD  modes  as  the  rank  of  POD  matrix  is  needed  to  reproduce  the  original 
CFD  solution; 

4.  The  eigenvalue  spectrum  of  the  ROM  stiffness  matrix  can  be  helpful  in  identifying  possible  numerically 
unstable  modes; 

5.  The  ROM  obtained  from  single -frequency  forcing  is  able  to  provide  reasonable  predictions  at  only  the 
forcing  frequency;  applying  forcing  functions  with  richer  frequency  content  provides  significant 
enhancement  in  the  predictive  capability  of  the  resulting  ROM’s. 


) 


■■K1[i+eT-0(x)]0’(x,T) 


III.  Formulation 


A.  Governing  equations 

The  governing  equations  are  the  quasi-one-dimensional  unsteady  Euler  equations  with  a  single-step  chemical 
reaction  and  a  specified  reaction  distribution: 


dQ  dE 

—  +  — 
dt  dx 


=  H  +  H/+H^ 


(2) 


where, 
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Here  x  and  t  are  the  space  and  time  variables,  p  is  the  density,  u  is  the  velocity,  e  is  the  total  energy,  p  is  the 
pressure,  Yox  is  the  oxidizer  mass  fraction  and  A  =  A(x)  is  the  cross-section  area  of  the  geometry.  The  effects  of  fuel 
addition  are  accounted  through  steady  source  term  H y,  where  (bf  =  Cf/0cbox  with  constant  Cf/0  representing  hiel-to- 
oxidizer  ratio  and  a  sinusoidal  spatial  distribution  is  used  to  model  the  oxidizer  reaction, 


co  -krpY 

ox  f  r  ox 


f  f  7  v* 

V  v  2  w.  „ 


1  +  sin 


,  (  V  ls  <  x  <  lf  ),  where  ls  and  lf  are  the  axial  locations  of  the  beginning 


and  end  of  the  combustion  zone.  The  reaction  constant  k  f  is  selected  to  ensure  that  the  oxidizer  is  consumed  within 
the  specified  combustion  zone.  The  unsteady  combustion  response  is  accounted  in  the  source  term,  Hq,  using  the  n-r 
,m  —  ,  .p(x,t-T)-p(x) 

model  [17],  q  =q-n-a[x) - - - - ,  which  relates  the  unsteady  heat  release  to  the  pressure  oscillations 

P\x) 

through  an  index,  n  and  a  time  lag  constant,  r.  Here  q  is  the  integral  mean  combustion  heat  release  and  a(x)  is  a 
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scaling  function  following  a  normal  distribution,  a(x)  =  — .^exp 


(t 

2cr2 


gA2k 

v  y 

[18]  to  simulate  combustion  instability  in  a  longitudinal  rocket  combustor, 

For  simplicity,  the  linearized  version  of  Eq.  (2)  is  used  for  the  studies  in  ROM  development, 


This  model  was  previously  used 
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Fluctuating  conditions  inside  the  computational  domain  are  obtained  by  specifying  a  periodic  upstream  mass 
flow  rate,  m’[t^ ,  oscillating  around  the  mean  value  at  the  inlet,  m0 , 


Nf 


m 


l(t)  =  m0—ZM24fo  +  V)t) 

y  k=l 


(4) 


where  /0  is  the  initial  frequency;  A/  represents  the  fundamental  frequency  increment  and  Nf  is  the  total  number  of 
frequencies  included  in  the  forcing  function.  The  fundamental  period  Tp  is  determined  by  A f  such  that  Tp  =1/A f. 
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Note  that  if  A/=/0,  Eq.  (4)  represents  a  standard  Fourier  series  although  herein  we  generally  take  A/ <  f0  so  that  Eq. 
(4)  differs  slightly  from  a  Fourier  series. 


B.  Construction  of  POD  eigen-bases  for  vector  equations 

POD  eigen-bases  are  calculated  based  on  the  CFD  solutions,  (x,t)  obtained  from  Eq.  (3), 


r  t,n  (X)  ) 

r  t,n  (X)  ' 

Np  Np 

t"(4q;  44  ~  (0cr»a>»  (x) = ZA 

n= 1  n= 1 
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4,44 

np 

=1X4 
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where  on  is  the  singular  value  of  the  nth  POD  mode  (scalar),  an(t)  is  the  nth  POD  temporal  mode  (scalar),  and  the  nth 
eigen-mode,  0„(x),  is  an  orthonormal  vector  function, 


1,  if  k  =  n 
0,  otherwise’ 


X  =  {0  <x<L] 


(6). 


It  should  be  noted  that,  unlike  in  the  scalar  equation  case,  to  obtain  reasonably  scaled  POD  eigen-bases  for  the 
vector  variables,  a  normalization  matrix  P(x)  must  be  used  before  calculating  the  POD  eigen-bases.  Likewise,  to 
reconstruct  the  CFD  solutions  using  the  POD  eigen-bases,  the  matrix  P(x)  again  needs  to  be  included, 


Q'p(x,t)*p-'  (4Zfl4) 


r  t,n  (4 A 

(4 

0T,n  (4 
KAx)  j 


(7). 


Definitions  of  the  normalization  matrix,  P(x),  are  considered  later. 


C.  Model  reduction  of  Euler  equations 

The  application  of  the  POD-Galerkin  method  to  the  linearized  Euler  equations  (Eq.  (3))  is  briefly  introduced  here. 
Additional  details  can  be  found  in  Ref  [11].  Upon  obtaining  the  eigen-bases  as  in  Eq.  (5),  the  target  governing 
equation  is  projected  onto  the  kth  eigen-mode  0^(x)  throughout  the  whole  computational  domain.  Before  the 
projection,  the  governing  equation  needs  to  be  normalized  by  pre-multiplying  by  the  matrix  P(x), 


(444 


eQ;  (x,t)  +  _ ^8Ap  (x)q;  (x,?4 


dt 


dx 


dx  =  \<s>l  (x)p(x)r;  {x)Dp  (x)Q;  (x,t)dx  (8), 


Substituting  the  POD  expansion,  Eq.  (7)  into  Eq.  (8)  and  using  a  numerical  quadrature  to  approximate  the  integrals, 
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(9). 


Following  the  model  reduction  procedure  in  Ref  [1 1]  and  using  a  consistent  discretization  scheme  to  approximate 
the  gradient  term  in  Eq.  (9)  as  discussed  in  Section  II,  an  ODE  system  can  be  obtained  with  the  contributions  from 
boundary  conditions  appearing  as  a  source  term  on  the  right-hand-side, 


M  -  £a(f )  =  h  (2) 

NI 

«*(°)=2>Lq  'p  (xz- , t  =  0)  Axt  (initial  conditions) 

2=1 


(10) 


where  a  (2)  =  [a,  (2)  •••  ak(t )  •••  aNf  (2)J  , 

h  (2)  =  [h}  (2)  •••  hk  (2)  •••  hN  (2)J  with  contributions  from  boundary  conditions,  hk  (2)  =  Fk  (w' (2))  , 

Z  is  the  stiffness  matrix  which  describes  the  dynamics  of  the  reduced  ODE  system. 


IV.  Overview  of  T est  Problem 

A  constant-area  duct  is  set  up  for  the  ROM  development  as  shown  in  Fig.  1  with  an  acoustically  open  boundary 
condition  at  the  outlet  and  a  mass  flow  perturbation  specified  at  the  inlet.  For  a  typical  computation,  the  steady  state 
numerical  solution  of  Eq.  (2)  is  computed  first  with  a  constant  mass  flow  rate,  rh0  at  the  inlet.  Representative  steady 
state  solutions  are  shown  in  Fig.  2(a),  the  temperature  of  which  has  the  general  character  of  an  Arrhenius  term  with  a 
rapid  increase  followed  by  an  asymptotic  approach  to  the  final  temperature.  The  upstream  boundary  condition  is 
then  switched  to  a  temporally  varying  periodic  condition  as  defined  in  Eq.  (4)  and  the  time-accurate  computation  is 
continued  until  stationary  conditions  are  reached. 

x/L  =  0  x/L  —  1 

m(t)  =  Mq  +m'(t)  =  m0 

s-  0.1,  /2back  -  0.45  MPa 


1+^-Zsin(M/o+A/)f) 


/^back 


Figure  1  Computational  setup  for  the  one-dimensional  Euler  equations. 
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The  solutions  are  tabulated  and  stored  periodically  (i.e.,  at  several  time-levels  within  a  forcing  period)  thereby 
generating  a  rectangular  matrix  that  can  be  used  as  a  database  for  fitting  eigen-bases  by  means  of  the  POD  procedure 
detailed  below.  The  POD  eigen-bases  are  then  applied  to  the  governing  linear  or  nonlinear  PDE  to  derive  the 
reduced-order  ODE  formulation  (or  ROM). 

Before  delving  into  the  POD  eigen-bases  generation  and  reduced-order  model  development,  it  is  important  to  get 
a  basic  understanding  of  the  inherent  system  response  of  the  test  problem.  For  this  purpose,  no  unsteady  combustion 
response  is  included  in  the  test  problem  (i.e.  n  is  set  to  0  in  Eq.  (3)).  To  identify  the  system  response,  several  sets  of 
CFD  solutions  are  calculated  using  different  mass  flow  forcing  frequencies  at  the  inlet  following  the  FTF/FDF 
(Flame  Transfer/Describing  Function)  approach  [19].  The  perturbation  level  s'm  Eq.  (4)  is  set  to  0.1.  For  each  case, 
only  one  frequency  is  included  in  the  forcing  function  and  the  amplitudes  and  phase  relations  with  regard  to  the 
boundary  perturbations  are  then  calculated  at  different  locations  accordingly.  The  resulting  amplitudes  and  phases  of 
pressure  (p’)  and  velocity  ( u  ’)  fluctuations  are  plotted  versus  the  forcing  frequencies  and  shown  in  Fig.  2(b).  A 
strong  resonant  system  response  at  approximately  1350Hz  can  be  clearly  seen,  which  corresponds  to  the  1st 
fundamental  acoustic  frequency  (IF)  determined  by  the  duct  length  and  mean  speed  of  sound.  We  point  out  that  this 
natural  system  response  is  a  unique  characteristic  of  the  Euler  equations  and  does  not  appear  in  the  scalar  equation 
case.  This  is  due  to  the  simultaneous  presence  of  both  forward  and  backward  moving  waves  in  the  Euler  system, 
which  naturally  gives  rise  to  resonant  (standing)  wave  solutions  that  correspond  to  the  fundamental  acoustic  modes. 


(a)  Steady  state  solutions  to  Eq.  (2) 


(b)  Generic  system  responses  (top:  limit  cycle 
amplitude;  bottom:  phase  angle) 


Figure  2  Characteristics  of  the  (a)  Steady  state  solutions  to  Eq.  (2)  and  (b)  Identified  generic  system 

responses  of  the  test  problem. 
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V.  Results 


A.  POD  Verification 


The  effectiveness  of  POD  eigen-bases  for  representing  the  original  numerical  solutions  is  verified  in  the  same 
way  as  in  the  scalar  equation.  At  the  same  time,  we  also  validate  the  two  bases  per  frequency  rule  identified  in  the 
scalar  equation  studies  (see  conclusion  1  in  Section  II).  The  cases  used  for  POD  verification  are  summarized  in 
Table  1.  Three  cases  are  selected  with  different  forcing  functions  from  single  to  multiple  frequencies  and  POD 
eigen-bases  are  generated  for  each  case.  As  discussed  earlier  in  regard  to  Eq.  (5),  the  CFD  solutions  to  the  Euler 
equations  need  pre-treatment  before  the  POD  eigen-basis  calculation.  In  this  section,  the  four  variables  ( u  \  p  \  T ’ 
and  Y’ox)  are  normalized  by  the  maximum  values  of  the  respective  fluctuating  quantity  throughout  the  whole  domain 
(i.e.  P(x)  =  diag (l  /  p^x ,  1  /  u'm!0l  ,1/T^,1/  Y'ox  m^  )  where  p'm  =  Max  {p'(x,t)} ,  V  0  <  *  <  L  and  t  etc.).  It  can  be 
readily  seen  that  the  two  bases  per  frequency  rule  is  still  valid  for  the  linearized  Euler  equations. 


Case  No. 

/o,Hz 

¥,  Hz 

#  of  forcing  frequencies 

Rank  of  POD  matrix 

G1 

1000 

1 

2 

G2 

1000 

250 

3 

6 

G3 

750 

250 

5 

10 

Table  1  Summary  of  cases  used  in  studies  of  POD  verifications. 


Similar  to  the  scalar  equation  studies,  the  test  solutions  are  reconstructed  using  the  POD  eigen-bases  as  given  in 
Eq.  (7)  and  then  compared  with  the  numerical  solutions  using  the  L2-norm  defined, 

J 1 4  f  >  0  -  pI  (*, )  Z  ■ an  (t ) (*, )  |  dt 

L2  -norm  (AM  of  variable  qk  at  x.  =  — - - t - - -  (11), 

\q'k1(xi,t)dt 

where  (xt )  is  the  kth  row  of  matrix  P~l  (x: ) .  The  L2-norms  are  plotted  for  each  case  in  Fig.  3  against  the  number 

of  POD  modes  included  in  the  solution  reconstruction  of  all  four  variables  ( u\p\  T’  and  7'OJC).  For  all  three  cases, 
the  L2-norms  saturate  and  reach  good  accuracy  after  including  a  mode  number  equal  to  the  matrix  rank  thereby 
validating  the  effectiveness  of  the  POD  eigen-bases  for  capturing  the  original  dynamics.  The  verified  POD  modes 
here  will  be  carried  on  next  to  build  the  ROMs  in  following  section.  So  far,  both  scalar  and  Euler  equation(s)  studies 
show  consistency  in  the  POD  eigen-bases  generation. 
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Figure  3  POD  verifications  of  cases  in  Table  1  (L2-norm  calculated  at  x!L  =  0.5). 


B.  ROM  Characteristics  for  Vector  Equations 

Knowing  the  number  of  POD  modes  needed  for  ROM  development  in  previous  section,  the  characteristics  of  the 
ROM  are  studied  in  this  section  for  the  three  cases  in  Table  1.  The  eigen- value  spectra  of  the  ROM  stiffness 
matrices  ( L  in  Eq.  (10))  are  investigated  and  shown  in  Fig.  4.  In  all  cases,  the  number  of  POD  modes  included  in  the 
ROM  development  is  equal  to  the  rank  of  POD  matrix.  For  Cases  G1  and  G3,  all  eigen-values  have  negative  real 
parts  (cr)  which,  as  discussed  in  the  scalar  equation  studies,  indicates  that  the  corresponding  ROM’s  will  be  stable. 
Solving  the  ROM  equation  for  these  two  cases  using  the  same  forcing  as  used  in  the  CFD  solution  gives 
reconstructed  solutions  that  are  indeed  in  excellent  agreement  with  the  CFD  results  as  shown  in  Fig.  5.  The  ROM’s 
are  able  to  reproduce  the  original  CFD  solutions  and  show  no  unphysical  or  unstable  behavior. 

Though  ROM  development  for  Case  G1  and  G3  follows  the  knowledge  learned  from  the  scalar  equation  studies, 
Case  G2  shows  some  differences.  An  unstable  mode  is  detected  in  Case  G2  (highlighted  by  the  red  circle  in  Fig.  4) 
even  though  all  necessary  POD  modes  are  included  to  build  the  ROM.  This  unstable  mode  leads  to  unphysical 
growth  in  the  ROM  solution  as  shown  in  Fig.  6.  The  growth  rate  of  the  unstable  eigenvalue  is  small,  enabling  the 
ROM  solution  to  compare  well  with  the  CFD  solution  at  early  times  (0.025  to  0.07s  shown  in  the  zoomed-in  box  in 
the  p  ’  plot)  but  the  ROM  prediction  goes  off  track  at  longer  times  as  the  unstable  mode  starts  to  dominate.  It  should 
be  noted  that  the  POD  modes  used  for  the  ROM  development  are  calculated  using  CFD  solutions  corresponding  to 
one  cycle  (defined  based  on  A/ in  Table  1),  the  same  as  in  the  scalar  case;  here,  CFD  solutions  from  0.025  to  0.029s 
are  selected  and  so  the  ROM  solution  is  technically  valid  based  on  the  information  given.  However,  the  appearance 
of  the  unstable  mode  strongly  impacts  the  ROM’s  predictive  capability. 
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T,  K  U,  m/s 


Figure  4  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  Cases  G1  to  G3  in  Table  1. 
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Figure  5  ROM  and  CFD  solutions  comparison  atx/Z  =  0.5  for  Cases  G1  to  G3  in  Table  1. 
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Figure  6  ROM  and  CFD  solutions  comparison  at  x!L  =  0.5  for  Cases  G2  in  Table  1. 


C.  Effects  of  normalization  methods  on  ROM  characteristics 

To  investigate  the  causes  of  the  unstable  mode  seen  in  Fig.  4,  three  different  normalization  methods  are  tried  for 
Case  G2  to  assess  their  effects  on  the  numerical  stability  of  the  ROM.  The  methods  used  for  these  studies  are 
summarized  in  Table  2.  Method  III  is  the  one  that  was  used  previously  for  the  results  shown  in  Figs.  4  -  6  above. 


Variables 

I 

II 

III 

u ’ 

X 

► 

A 

P’ 

X 

► 

A 

T 

X 

► 

A 

Y’ 

J  ox 

X 

X 

A 

x:  no  normalization;  ► :  domain  averaged  mean  value,  e.g.  pKi  = 


l-\Lap(x)dx 


A:  Maximum  value  of  all  points,  e.g.  pref  =  Max  { p'(x,t)},  V  0  <x<L  and  t 


Table  2  Summary  of  normalization  methods  used  in  studies  of  ROM  stability  issues  for  Case  G2. 
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(b)  Method  II 


(c)  Method  III 

Figure  7  Comparisons  of  normalized  CFD  solutions  using  methods  in  Table  2. 


Before  proceeding  to  investigations  of  normalization  effects  on  the  ROM  characteristics,  the  normalized 
solutions  used  for  POD  eigen-bases  generation  are  compared  for  the  different  methods  in  Fig.  7.  As  expected,  with 
no  normalization  (Method  I),  the  dominance  of  p  ’  overwhelms  the  other  three  variables  ( u  ’,  T’  and  Y’ox).  By  using 
mean  values  (Method  II),  amplitudes  of p  \  T’  and  Y ’ox  are  comparable  while  u  ’  contributes  more  than  the  other  three 
variables  in  POD  eigen-bases  generation  by  a  factor  of  approximately  10  so  the  mode  shapes  of  resulting  POD 
eigen- functions  are  expected  to  be  governed  by  the  velocity  fluctuations.  Method  III  seems  to  be  optimum  in 
keeping  four  variables  within  similar  fluctuating  levels,  which  makes  the  contributions  from  each  variable  to  the 
POD  eigen-bases  generation  comparable. 
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Figure  8  Comparisons  of  ROM  eigen-values  using  different  normalization  methods  in  Table  2. 


An  overview  of  the  developed  ROMs’  eigen- values  for  the  different  normalization  methods  is  given  in  Fig.  8 
with  unstable  modes  highlighted  by  red  circles.  First,  it  can  be  seen  that  the  ROM  can  be  stabilized  by  using  Method 
II,  which  enables  the  ROM  to  reproduce  the  original  CFD  solutions  as  shown  in  Fig.  9  and  no  unphysical  or  unstable 
behavior  seen  in  Fig.  6  is  present.  Second,  with  no  normalization,  the  resulting  ROM  shows  a  very  large  cr  value, 
which  can  lead  to  unreasonable  growth  and  totally  derail  the  ROM  predictions. 
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Figure  9  ROM  and  CFD  solutions  comparison  at  x!L  =  0.5  for  Cases  G2  using  Method  II  in  Table  2. 
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Having  determined  that  the  unstable  mode  identified  in  Fig.  4  can  be  avoided  by  using  a  different  normalization 
method  for  Case  G2,  we  next  test  the  validity  of  Method  II  for  the  other  two  cases  in  Table  1.  The  corresponding 
ROM  eigen- value  spectrum  is  shown  in  Fig.  10.  It  can  be  seen  that  the  ROM  for  Case  G1  is  still  stable  using 
Method  II  for  normalization.  However,  an  unstable  mode  (red  circled)  can  be  observed  in  the  ROM  for  Case  G3  at 
the  highest  frequency  (~3200Hz)  with  a  large  cr  value  (>600s_1).  Based  on  the  observation  here,  it  seems  that  a  valid 
normalization  method  for  one  case  does  not  necessarily  guarantee  a  universal  numerical  stability  for  the  other  cases. 
Moreover,  the  effects  of  normalization  methods  seem  to  be  case-dependent  and,  in  fact,  are  reliant  on  the  forcing 
function,  which  determines  the  spatial  mode  shapes  of  CFD  solution.  In  is  noteworthy  that  Method  II  fails  for  Case 
G3,  which  only  involves  two  additional  frequencies  (750  and  1750Hz)  compared  to  the  stable  Case  G2. 


Figure  10  ROM  eigen-value  spectrum  for  Cases  G1  and  G3  using  Method  II. 


D.  Effects  of  Normalization  Methods  on  ROM  Predictive  Capability 

After  obtaining  some  basic  understanding  of  how  normalization  methods  can  affect  ROM  stability,  the  predictive 
capabilities  of  ROMs’  developed  using  different  normalization  methods  are  assessed  in  this  section.  Studies  are 
performed  for  ROMs  developed  using  both  single-  and  multi-frequency  forcing. 


Single-frequency  forcing 

Two  single-frequency  forcing  functions  are  used  to  train  ROMs  for  predictive  capability  investigations  here,  one 
of  which  is  at  1000Hz  (Case  G1  in  Table  1)  and  the  other  one  is  chosen  to  be  1350Hz  (denoted  as  Case  G4  here), 
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which  corresponds  to  the  natural  system  frequency  identified  in  Fig.  2(b).  The  studies  here  are  focused  on  how 
variable  normalization  affects  the  resulting  ROM’s  predictive  capability. 

For  both  Case  G1  and  G4,  Method  II  and  III  are  chosen  for  ROM  predictive  capability  evaluations.  The  ROM 
predictive  capability  is  evaluated  in  a  different  way  compared  to  the  previous  scalar  equation  studies.  Here,  the 
ROMs  are  used  to  predict  the  limit  cycle  amplitude  and  phase  angle  at  different  input  frequencies,  which  are  then 
compared  with  CFD  results  as  shown  in  Fig.  10.  Consistent  with  observations  in  the  scalar  equation  studies,  the  two 
ROMs  are  able  to  predict  the  correct  amplitude  and  phase  angle  at  the  forcing  frequencies  for  both  cases.  The  Case 
G1  ROMs  are  only  able  to  capture  responses  at  1000Hz  while  predictions  at  other  frequencies  are  poor,  which  is 
expected  from  scalar  equation  studies  (Conclusion  5  in  section  II).  Applying  different  normalization  methods  for 
Case  G1  does  not  change  the  overall  ROM  performance. 

The  Case  G4  ROMs  show  more  capabilities  for  predicting  frequency  response  not  included  in  the  forcing 
function.  The  Method  II  ROM  is  able  to  capture  the  overall  trend  of  the  system  response,  but  seems  to  predict  a 
faster  amplitude  decay  versus  frequency  than  the  CFD  solution  does.  Interestingly,  switching  the  normalization 
method  seems  to  have  more  significant  effects  than  in  Case  Gl.  The  Method  III  ROM  shows  a  sharper  resonant 
peak  than  Method  II  in  limit  cycle  amplitude,  which  can  also  be  reflected  by  an  even  faster  amplitude  decay  at 
frequencies  not  equal  to  1350Hz.  This  might  be  attributed  to  the  fact  that  Method  III  puts  more  comparable 
weightings  on  each  variable  in  generating  POD  eigen-bases  while,  in  Method  II,  the  eigen-bases  for  ROM 
development  are  dominated  by  velocity  fluctuations  ( u  ’),  but  further  investigations  are  needed. 


(a)  Case  Gl  ROM  (1000Hz  forcing)  (b)  Case  G4  ROM  (1350Hz  forcing) 


Figure  11  Comparisons  of  ROM  predictive  capability  using  different  normalization  methods  at  x/L  = 

0.5  for  Case  Gl  and  G4. 
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Multi-frequency  forcing 

The  next  set  of  evaluations  is  carried  out  for  ROMs  trained  using  multi-frequency  boundary  forcing  (Case  G2  & 
G3  in  Table  1).  Method  II  and  III  are  used  for  comparisons  to  assess  the  effects  of  normalization  strategies  for  the 
vector  system,  which  are  shown  in  Fig.  11.  It  has  already  been  shown  in  Fig.  8  that  Method  II  generates  stable  ROM 
while  Method  III  can  cause  numerical  instability  issues  for  Case  G2,  which  leads  to  a  total  failure  in  the  ROM 
predictions  of  limit  cycle  amplitude  and  phase  angle,  while  the  stable  ROM  (from  Method  II)  is  able  to  capture  the 
overall  CFD  solution  trend  as  shown  in  Fig.  12(a). 

For  Case  G3,  on  the  other  hand,  we  have  noted  that  Method  II  is  unstable  and  Method  III  is  stable  and  the  results 
here  confirm  that  Method  III  is  able  to  obtain  very  good  agreement  with  the  CFD  solutions.  Further,  by  comparing 
the  valid  ROM  from  each  case  (Method  II  for  Case  G2  and  Method  III  for  Case  G3),  it  can  be  seen  that  by  adding 
two  more  frequencies  (750  and  1750Hz)  in  the  forcing  function,  the  resulting  ROM  predictive  capability  can  be 
improved  a  lot,  which  is  consistent  with  Conclusion  5  from  scalar  equation  studies,  and  confirms  the  value  of 
applying  forcing  functions  with  richer  frequency  content  in  the  ROM  development. 


(a)  Case  G2 


(b)  Case  G3 


Figure  12  Comparisons  of  ROM  predictive  capability  using  different  normalization  methods  for  Case  G2 

&  G3  in  Table  1  at  x!L  =  0.5. 
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VI.  Conclusions 


Explorations  of  the  POD/Galerkin  method  for  developing  a  reduced  order  model  (ROM)  for  the  linearized 
reactive  Euler  equations  are  performed  as  an  extension  of  previous  studies  using  a  reaction-advection  scalar 
equation.  Consistent  with  the  previous  scalar  equation  study,  the  approach  utilizes  CFD  solutions  to  obtain  the  POD 
basis  vectors,  which  are  then  used  within  a  Galerkin  formulation  to  reduce  the  governing  partial  differential 
equations  to  a  set  of  ordinary  differential  equations  (ODE).  The  studies  in  this  paper  are  focused  on  two  objectives: 
first,  to  validate  the  conclusions  from  the  scalar  equation  study  for  the  Euler  equations;  second,  to  identify  specific 
issues  present  in  ROM  development  for  vector  equations  as  a  background  for  future  investigations  using  the  full 
Navier- Stoke  equations. 

In  Euler  equations  studies,  consistent  conclusions  can  be  reached  as  in  the  scalar  equation  case  with  regard  to  the 
number  of  POD  modes  required  for  ROM  construction  and  the  enhancement  of  the  predictive  capability 
enhancement  by  enriching  the  frequency  content  of  the  forcing  function  used  to  train  the  ROM.  In  other  words,  we 
are  able  to  confirm  that  the  two  bases  per  frequency  rule  remains  operative  for  the  Euler  equations  and  the  addition 
of  multiple  frequencies  improves  the  ROM  prediction  for  a  broader  frequency  band. 

The  major  issue  identified  in  the  ROM  development  of  vector  equations  is  the  need  to  normalize  the  vector 
variables  used  to  generate  the  coupled  POD  eigen-bases,  which  is  not  needed  in  the  scalar  equation.  Without 
normalization,  we  observe  strong  unstable  behavior  of  the  ROM  eigenvalue  spectrum,  rendering  the  corresponding 
ROM  incapable  of  predicting  the  correct  system  response.  Normalizing  the  variables  by  the  maximum  value  in  the 
domain  provides  stable  behavior  for  some  cases,  but  unstable  modes  are  observed  for  other  cases.  Alternate 
normalizations  such  as  using  the  mean  value  of  the  variables  stabilizes  some  of  the  previously  unstable  cases,  but 
destabilizes  other  stable  cases.  The  unstable  modes  are  shown  to  lead  to  numerical  stability  problems,  which  are 
characterized  by  unphysical  amplitude  growth  in  the  ROM  predictions,  while  the  correct  CFD  solution  reaches  a 
limit  cycle  (i.e.,  no  growth). 

Overall,  applying  ROM  capabilities  to  the  prediction  of  combustion  instability  remains  an  important,  but  highly 
challenging  task.  Such  models  have  the  potential  to  provide  effective  predictions  at  a  fraction  of  the  cost  needed  for 
solving  the  original  partial  differential  equation  set.  However,  as  the  present  paper  points  out,  significant  issues 
remain  in  formulating  the  ROM  for  vector  systems  such  as  the  Euler  equations.  More  studies  need  to  be  done  to 
understand  the  effects  of  variable  normalization  on  the  ROM  construction  and  the  proper  normalization  procedure  to 
avoid  spurious  unstable  modes.  A  possible  solution  is  to  re-defme  the  inner  product  needed  for  the  coupled  POD 
eigen-bases  calculation  so  that  the  coupling  between  different  unsteady  quantities  can  be  incorporated  naturally  in 
POD  modes  generation.  An  example  of  that  can  be  found  in  Ref.  [16],  which  will  be  a  starting  point  for  our  future 
studies. 
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